##########################################################################################

library(Seurat)
library(data.table)
library(optparse)
library(ggplot2)
library(dplyr)
library(RColorBrewer)

##########################################################################################

option_list <- list(
    make_option(c("--singleCell_sample_file"), type = "character"),
    make_option(c("--single_cell_file"), type = "character"),
    make_option(c("--out_path"), type = "character")
)

if(1!=1){

    singleCell_sample_file <- "~/20220915_gastric_multiple/dna_combinePublic/config/singleCell_Sample.useThree.list"
    single_cell_file <- "~/20220915_gastric_multiple/dna_combinePublic/public_ref/singleCell/njmu/ALL_SIN_celltype.Rdata"
    out_path <- "~/20220915_gastric_multiple/dna_combinePublic/finalPlot/revise/im_favored"

}

##########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

singleCell_sample_file <- opt$singleCell_sample_file
single_cell_file <- opt$single_cell_file
out_path <- opt$out_path

dir.create(out_path , recursive = T)

##########################################################################################

info_singlecell <- data.frame(fread(singleCell_sample_file))
sc_dataset <- load(single_cell_file, verbose = F)

##########################################################################################

sc_dataset_all <- ALL_SIN_celltype

##Idents函数定义要取的类别
Idents(sc_dataset_all) <- "sample"   

##########################################################################################
## 提取最后纳入分析的样本
## 3个
info_singlecell <- subset( info_singlecell , singlecell_ID != "" )
patientid <- unique(substring(unique(info_singlecell$singlecell_ID) , 0 , 5))
## 提前用到样本的细胞
im_sc_dataset <- subset(sc_dataset_all , patient %in% patientid & sample == 'IM' & celltype %in% c("Bcell" , "Mast" , "Tcell"))
im_sc_dataset$MUC6_mut <- ifelse( im_sc_dataset$patient == "JZ732" , "Mut" , "Wild" )

mutTumor <- "JZ732"
wildTumor <- patientid[!patientid %in% mutTumor]


##########################################################################################

trans <- function(num){
    up <- floor(log10(num))
    down <- round(num*10^(-up),2)
    text <- paste("P == ",down," %*% 10","^",up)
    return(text)
}

##########################################################################################
## 按照不同病理类型提取
## 柱状图
## 合并在一起
sc_dataset <- im_sc_dataset

## 野生型样本
wild_cells <- sc_dataset$celltype1[grep( paste0(wildTumor , collapse = "|") , names(sc_dataset$celltype) )]

## 突变型样本
mut_cells <- sc_dataset$celltype1[grep( paste0(mutTumor , collapse = "|") , names(sc_dataset$celltype) )]

## 统计细胞个数
wild_cells <- data.frame(table(wild_cells))
mut_cells <- data.frame(table(mut_cells))


##########################################################################################
## 统计免疫细胞
# CD1C: CD1C是CD1分子的一种，CD1分子是一类抗原呈递分子，它们可以呈递脂质抗原给T细胞，从而参与免疫监视。
# CD8: CD8是T细胞的一个亚群，这些细胞具有细胞毒性，能够识别并杀死被病毒感染的细胞或癌细胞，是免疫监视的关键执行者。
# Mast: 肥大细胞（Mast cells）在免疫监视中发挥作用，尤其是在对寄生虫感染的早期反应中，它们可以释放炎症介质。
# M1: M1型巨噬细胞是一种经典活化的巨噬细胞，它们具有强烈的抗病原体和抗肿瘤活性，参与免疫监视。
# M2: M2型巨噬细胞通常与组织修复和炎症调节有关，但它们也可以在免疫监视中发挥作用，尤其是在清除病原体和促进肿瘤生长的环境中。
# Plasma: 浆细胞是B细胞的终末分化形式，它们产生并分泌抗体，这些抗体可以识别并结合到病原体上，从而标记它们供其他免疫细胞消灭。
# Th17: Th17细胞是一种辅助性T细胞，它们在对抗细菌和真菌感染中发挥重要作用，尤其是在粘膜表面。
# Tcm: T细胞中央记忆细胞（Tcm）是T细胞的一种，它们在免疫监视中起到维持长期免疫记忆的作用。
# Treg: 调节性T细胞（Treg）在免疫监视中起到调节免疫反应的作用，防止免疫过度激活，但也可能有助于肿瘤逃避免疫监视。
# Naive幼稚T细胞是那些尚未遇到其特异性抗原的T细胞。它们是T细胞发育过程中的一个早期阶段，尚未被激活以对抗特定的病原体

wild_cells$Ratio <- wild_cells$Freq/sum(wild_cells$Freq)
mut_cells$Ratio <- mut_cells$Freq/sum(mut_cells$Freq)

wild_cells$Type <- "Wild"
mut_cells$Type <- "Mut"
colnames(wild_cells)[1] <- "Cell_Type"
colnames(mut_cells)[1] <- "Cell_Type"

dat_plot <- rbind(wild_cells , mut_cells )

result_plot <- c()
for( cell_Type in unique(dat_plot$Cell_Type) ){
        dat_tmp <- subset( dat_plot , Cell_Type == cell_Type )
        a <- round(100 * subset( dat_plot , Cell_Type == cell_Type & Type == "Mut"  )$Ratio)
        if(length(a)==0){a=0}
        b <- 100 - a

        c <- round(100 * subset( dat_plot , Cell_Type == cell_Type & Type == "Wild"  )$Ratio)
        if(length(c)==0){c=0}
        d <- 100 - c

        p <- fisher.test( matrix(c(a,b,c,d) , ncol = 2) )$p.value

        dat_tmp <- data.frame( 
                Cell_Type = c("Cells of interest" , "Other Cells" , "Cells of interest" , "Other Cells") ,
                Ratio = c(
                        subset( dat_plot , Cell_Type == cell_Type & Type == "Mut"  )$Ratio ,
                        1 - subset( dat_plot , Cell_Type == cell_Type & Type == "Mut"  )$Ratio ,
                        subset( dat_plot , Cell_Type == cell_Type & Type == "Wild"  )$Ratio ,
                        1 - subset( dat_plot , Cell_Type == cell_Type & Type == "Wild"  )$Ratio
                        ),
                Type = c("Mut" , "Mut" , "Wild" , "Wild") ,
                Cell_Type_use = cell_Type
         )
        if( p < 0.01 ){
                p_text <- trans(p)
        }else{
                p_text <- paste0( "P == " , round(as.numeric(p) , 3) ) 
        }
        dat_tmp$p <- ""
        dat_tmp$p_text <- ""
        dat_tmp$p[1] <- p
        dat_tmp$p_text[1] <- p_text
        result_plot <- rbind( result_plot , dat_tmp )
}

result_plot$value_text <- paste0( round(result_plot$Ratio , 2) * 100 , "%") 
result_plot$value_text <- ifelse( round(result_plot$Ratio , 2) == 0 , "" , result_plot$value_text  )

col <- c(
    rgb(red=179,green=34,blue=35,alpha=255,max=255), 
    rgb(red=2,green=100,blue=190,alpha=255,max=255) 
)
col <- col[2:1]
names(col) <- c("Other Cells" , "Cells of interest")

cell_order <- c("Plasma" , "Naive" , "CD8" , "Tcm" , "Th17" , "Treg" , "Mast")
result_plot <- data.frame(result_plot)
result_plot$Cell_Type_use <- factor( result_plot$Cell_Type_use , levels = cell_order , order = T )

p1 <- ggplot(result_plot,aes(x=Type,y=Ratio,fill=factor(Cell_Type))) +
        geom_bar(stat="identity") +
        ylab("Proportion") +
        geom_text(aes(label=p_text , y = 1.05 , x = 1.5),parse = TRUE,size=3.5 , color = "black") +
        geom_text(aes(label=value_text) , position=position_stack(vjust = 0.5) , size=4 , color="white")+
        xlab(NULL) +
        facet_grid( .~Cell_Type_use , scales = "free_x") +
        theme_bw() +
        theme(
                panel.grid.major=element_blank(),
                panel.grid.minor=element_blank(),
                panel.background = element_blank(),
                panel.border = element_blank(),
                legend.position ='top',
                legend.title = element_blank() ,
                legend.text = element_text(size = 12,color="black",face='bold'),
                axis.text.x = element_text(size = 12,color="black",face='bold'),
                axis.text.y = element_text(size = 10,color="black",face='bold'),
                axis.title.x = element_text(size = 10,color="black",face='bold'),
                axis.title.y = element_text(size = 14,color="black",face='bold'),
                strip.text.x = element_text(size = 12,color="black",face='bold'),
                axis.ticks.length = unit(0.2, "cm") ,
                axis.line = element_line(size = 0.5))  +
        scale_fill_manual(values=c(col))


out_name <- paste0( out_path , "/STAD_MUC6_mutwild.IM.Immune.bar.CellRate.pdf"  )
ggsave(file=out_name,plot=p1,width=10,height=4)

out_name <- paste0( out_path , "/STAD_MUC6_mutwild.IM.Immune.bar.CellRate.tsv"  )
write.table( result_plot , out_name , row.names = F , quote = F , sep = "\t" )

##########################################################################################
## 盒图展示每个样本
result_final <- c()

## 野生型样本
dat_wild_cells <- c()
for( tumor in wildTumor ){
        mut_cells <- im_sc_dataset$celltype1[ grep( tumor , names(im_sc_dataset$celltype)) ]
        if(length(mut_cells) > 0){
                mut_cells <- data.frame(table(mut_cells))
                mut_cells$Ratio <- mut_cells$Freq/sum(mut_cells$Freq)
                colnames(mut_cells)[1] <- "Cell_Type"
                mut_cells$Sample <- tumor
                dat_wild_cells <- rbind( dat_wild_cells , mut_cells )
        }
}
dat_wild_cells$Type <- "Wild"

## 突变型样本
dat_mut_cells <- c()
for( tumor in mutTumor ){
        mut_cells <- im_sc_dataset$celltype1[ grep( tumor , names(im_sc_dataset$celltype)) ]
        if(length(mut_cells) > 0){
                mut_cells <- data.frame(table(mut_cells))
                mut_cells$Ratio <- mut_cells$Freq/sum(mut_cells$Freq)
                colnames(mut_cells)[1] <- "Cell_Type"
                mut_cells$Sample <- tumor
                dat_mut_cells <- rbind( dat_mut_cells , mut_cells )
        }
}

if(length(dat_mut_cells) > 0){
        dat_mut_cells$Type <- "Mut"
        dat_plot <- rbind(dat_wild_cells , dat_mut_cells )

        result_plot <- c()
        for( cell_Type in unique(dat_plot$Cell_Type) ){
                dat_tmp <- subset( dat_plot , Cell_Type == cell_Type )

                a <- dat_tmp[dat_tmp$Type=="Mut","Ratio"]
                if(length(a)==0){a=0}
                b <- dat_tmp[dat_tmp$Type=="Wild","Ratio"]

                p <- wilcox.test( a , b )$p.value
                if( p < 0.01 ){
                        p_text <- trans(p)
                }else{
                        p_text <- paste0( "p == " , round(as.numeric(p) , 3) ) 
                }
                dat_tmp$p <- ""
                dat_tmp$p_text <- ""
                dat_tmp$p[1] <- p
                dat_tmp$p_text[1] <- p_text
                result_plot <- rbind( result_plot , dat_tmp )
        }

        result_final <- result_plot

        col <- c(
            rgb(red=179,green=34,blue=35,alpha=255,max=255), 
            rgb(red=2,green=100,blue=190,alpha=255,max=255) 
        )

        names(col) <- c("Mut" , "Wild")
        result_plot$Cell_Type <- factor( result_plot$Cell_Type , levels = cell_order , order = T )
        result_plot <- subset( result_plot , Cell_Type != "Tumor" )
        result_plot$size_dot <- ifelse( result_plot$Type == "Wild" , 1.2 , 1.3 )
        result_plot$Type <- factor( result_plot$Type , levels = c("Wild" , "Mut") , order = T )

        p1 <- ggplot(result_plot, aes(x = Cell_Type , y = Ratio )) +
        geom_boxplot(size = 1.2 , outlier.alpha=0 , color = col['Wild']) + ## 去除散点，加粗线
        geom_jitter(position = position_jitter(0.2) , aes(color = Type , shape = Type )  , size = 3) + 
        scale_color_manual(values=c(col)) +
        #facet_grid( .~Cell_Type , scales = "free_x") +
        #geom_text(aes(label=p_text , y = 0.75 , x = 1.5),parse = TRUE,size=4 , color = "black") +
        xlab(NULL) +
        ylab("Proportion")+
        theme_bw() +
        theme(
            legend.position = 'top',
            legend.title = element_blank() ,
            panel.grid.major=element_blank(),
            panel.grid.minor=element_blank(),
            panel.background = element_blank(),
            panel.border = element_blank(),
            plot.title = element_text(size = 8,color="black",face='bold'),
            legend.text = element_text(size = 8,color="black",face='bold'),
            axis.text.y = element_text(size = 8,color="black",face='bold'),
            axis.title.x = element_text(size = 8,color="black",face='bold'),
            axis.title.y = element_text(size = 8,color="black",face='bold'),
            axis.text.x = element_text(size = 8,color="black",face='bold') ,
            axis.ticks.length = unit(0.2, "cm") ,
            axis.line = element_line(size = 0.5)) 

        out_name <- paste0( out_path , "/STAD_MUC6_mutwild.IM.Immune.box.CellRate.pdf"  )
        ggsave(file=out_name,plot=p1,width=4,height=4)
}

out_name <- paste0( out_path , "/STAD_MUC6_mutwild.IM.Immune.box.CellRate.tsv"  )
write.table( result_final , out_name , row.names = F , quote = F , sep = "\t" )
